Thermodynamics of the two-component Fermi gas with unequal masses at unitarity 
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We consider mass-imbalanced two-component Fermi gases for which the unequal-mass atoms 
interact via a zero-range model potential with a diverging s-wave scattering length a„, i.e., with 
1 /a s = 0. The high temperature thermodynamics of the harmonically trapped and homogeneous 
systems are examined using a virial expansion approach up to third order in the fugacity. We find 
that the universal part of the third-order virial coefficient associated with two light atoms and one 
heavy atom is negative, while that associated with two heavy and one light atom changes sign from 
negative to positive as the mass ratio k increases, and diverges when Efimov physics sets in at 
k = 13.61. By examining the Helmholtz free energy, we find that the equilibrium polarization of 
the trapped and homogeneous systems is for k = 1, but finite for k / 1 (with a majority of heavy 
particles). Compared to the equilibrium polarization of the non-interacting system, the equilibrium 
polarization at unitarity is increased for the trapped system and decreased for the homogeneous 
system. We find that unequal-mass Fermi gases are stable for all polarizations. 

PACS numbers: 



I. INTRODUCTION 

It has been predicted that ultracold two-component 
Fermi gases with mismatching Fermi surfaces sup- 
port novel phases such as the Fulde-Ferrell-Larkin- 
Ovchinnikov PHI] state or other interesting phase sep- 
arations, both in the homogeneous and inhomogeneous 
systems 043- I n equal-mass systems, the two Fermi sur- 
faces differ if the number of fermions of species 1 and 2 
differ, i.e., if the system exhibits a spin imbalance [Tolfllj ]. 
Alternatively, mismatching Fermi surfaces arise if species 
1 and 2 have different masses [l2T - fl9l ]. Experimentally, 
mass-imbalanced systems can be realized by simultane- 
ously trapping, e.g., 6 Li and 40 K The appli- 
cation of an external magnetic field in the vicinity of 
a Fano-Feshbach resonance allows for the realization of 
strongly-interacting gases, thus motivating our studies of 
the unequal-mass Fermi gas at unitarity. 

At the few-body level, two-component unequal mass 
systems with short-range interspecies interactions and 
mass ratio n are interesting because they support a vari- 
ety of intriguing states [27[ ■ At unitarity, the three-body 
system consisting of two heavy and one light particle [the 
(m,n2) = (2,1) system] with zero-range interspecies in- 
teractions supports no three-body bound state in free 
space for k < 8.62. For 8.62 < n < 13.61, however, 
the (774, n 2 ) = (2,1) system can support a bound state 
whose binding energy is determined by a three-body pa- 
rameter 28-32]. For yet larger mass ratios (n > 13.61), 
the three-body system in free space supports an infinite 
number of three-body Efimov states, which are geomet- 
rically spaced [3314351 ] . The spacing depends on the mass 
ratio k between the heavy and light particles with the 
energy of the most deeply bound state being determined 
by the so-called three-body Efimov parameter [36[ . 

The objective of this paper is to investigate the ther- 
modynamics of inhomogeneous and homogeneous two- 
component Fermi gases as a function of the mass ratio k. 



Throughout, we consider the unitary regime where the 
interspecies s-wave scattering length diverges and deter- 
mine the system properties as functions of the polariza- 
tion P, P = (Ai - N 2 )/(N 1 + N 2 ), and the tempera- 
ture T. We are limiting ourselves to the so-called high 
temperature regime, where the virial equation of state 
up to the third order is expected to provide a valid de- 
scription [13, HS|. For the trapped equal mass system, 
the virial equation of state up to the third order has 
been shown to be applicable over the temperature range 
of well above the Fermi temperature Tp down to about 
T F /2 [Mil. 

Our key findings are: (i) The second order virial equa- 
tion of state in a harmonic trap is independent of the 
mass ratio n. The mass ratio first enters at the third or- 
der in the virial expansion, (ii) The equilibrium state of 
the trapped system at unitarity is stable and favors a ma- 
jority of heavy particles for k 7^ 1. (iii) The equation of 
state of the homogeneous system depends strongly on the 
mass ratio. The mass ratio first enters at the zeroth or- 
der in the virial expansion, (iv) The equilibrium state of 
the homogeneous system at unitarity is stable and favors 
a majority of heavy particles. The equilibrium polariza- 
tion exhibits an intricate dependence on the third-order 
virial coefficients. 

The paper is organized as follows. Section [TTJ defines 
the virial expansion and the first few virial expansion co- 
efficients for a two-component Fermi gas in a harmonic 
trap where both species feel the same angular trapping 
frequency u>. The third order virial expansion coefficients 
at unitarity are determined as functions of the temper- 
ature and the mass ratio. Sections Mil and IIVI explore 
the virial equation of state of the harmonically trapped 
and homogeneous two-component Fermi gas at unitar- 
ity. Finally, Sec. [V] concludes. Appendix [A] summarizes 
the virial expansion for the harmonically trapped single- 
component Fermi gas and connects the results with those 
for the two-component Fermi gas. 
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II. DETERMINATION OF THE VIRIAL 
COEFFICIENTS 

The model Hamiltonian H (711,712) that describes un- 
equal mass two-component Fermi gases with 71 1 heavy 
and 7i2 light particles in a harmonic trap with angular 
trapping frequency u) reads 



exp[fii/(k B T)], the thermodynamic potential reads 



00 00 



.ni—0 no— 



(•5) 



Throughout, we are interested in the large T limit 
where Zi is small. Taylor expanding Eq. (J5J) about zi — 
to third order, we find 41] 



(1) 



m+112 / ^2 



j=rii+l 



ni ni+ri2 



+ E f^ + r^J+E E ^«). 



where 



i=i j=m+i 



(i) 



where mi denotes the mass of the heavy species and 7712 
denotes the mass of the light species. We define the mass 
ratio /c as k — 7711/7772 (k > 1). In Eq. ([1}, 7^ denotes 
the position vector of the j th atom measured with re- 
spect to the trap center. The interspecies interaction 
potential Vtb is parameterized by a zero-range potential, 
V t b(nj) = {2irh 2 a s / 'nred)S(rij)(d/drij)rij, where fi red de- 
notes the reduced mass, \i re d = 77117712/ (rn\ +7712), a s the 
interspecies s-wave scattering length, and r,j the inter- 
particle distance, r.y = \r% — 

If the complete energy spectrum of the Hamiltonian 
given in Eq. ([T]) were known, one could calculate ther- 
modynamic quantities like the average energy U or the 
entropy S from the free energy F, F — — k B T In Q Tlli „ 2 , 
using the canonical ensemble. Here, k B is Boltzmann's 
constant and Q ni „ 2 is the canonical partition func- 
tion ss, 



III — 1 

OO 

= -k B TQ h0 b , 



! ,"2 ^2 



(6) 

(7) 
(8) 



ri2 = l 



and 



AQ = - k B TQ lj0 [ 6 M ziz 2 - 
61,2 Z\zl + 62,1 zfz 2 



(9) 



? ni ,n 2 = Tr exp[-H(n 1 ,n 2 )/(k B T)], (2) s i ng l e heavy particle of species 1. The b nu0 (*>o,i 



Equations ©-© show that the thermodynamic poten- 
tial of the two-component Fermi gas can be written as 
a sum of the thermodynamic potentials Cl^ and 
of components 1 and 2 and a correction term AQ. The 
thermodynamic potentials f2^ and fij 01 the single- 
component species 1 and 2 are written in terms of the 
expansion or virial coefficients b nit o and 60, n2 (see Ap- 
pendix [X| and the canonical partition function Qi,o of a 

with 



where Tr is the trace operator. To evaluate the trace, we 
insert a complete set of eigenstates, yielding 



(3) 



The summation index j collectively denotes the complete 
set of quantum numbers allowed by symmetry. 

In practice, it is easier to work in the grand canonical 
ensemble, which can be considered to be a collection of 
canonical ensembles with rii and 712 particles of species 
1 and 2 in thermal equilibrium with each other. In the 
grand canonical ensemble, the chemical potentials /ii and 
fi2 of the two species are fixed. This implies that the 
system is characterized by the average numbers Ni and 
N2 of atoms of species 1 and 2. The thermodynamic 
potential fl^ of the two-component Fermi gas in the 
grand canonical ensemble is [45[ 



ft (2) = - fc B TlnTr exp[- (#(711,712) 
- /iini - fJ,2n 2 )/(k B T)], 



(4) 



where the trace operator now extends over tii and 
712. Rewritten in terms of the fugacities z,-, z% = 



ni > 1 (712 > 1) arise from the fermionic nature of the 
atoms and can be viewed as corrections to the Boltzmann 
gas [44, 45] . The correction term Afi contains the expan- 
sion or virial coefficients b ni ^ n2 , which account for the 
interactions between distinguishable fermions [ill |47j , 

61,1 =(Qi,i - <3i,oQo,i)/Qi,o, (10) 

h,2 =(Qi,2 - <3i,oQo,2 - &i,i<2i,o<3o,i)/Qi,0) (n) 

and 

&2,i =(Q2,i - Q2.0Q0.1 - &i,iQi,o<9i,o)/Qi,o- (12) 

In the non-interacting limit, i.e., for a s — 0, we label 
the partition functions by a superscript NI. In this case, 
Qn\,n 2 reduces to Qm,oQo,n 2 , since H(ni,n 2 ) is separa- 
ble for a s = 0. In this limit, 61.1 vanishes since 
equals Qi,oQo,i- Consequently, &i : 2 and 62.1 as well as 
all higher order virial coefficients fe nii „ 2 with 71 1 + 71 2 > 3 
vanish; thus, Af2 is zero for vanishing a s . 

To determine the b nitTl2 , one must know three things: 

(i) all bij where i < 71 1 and j < 712 or i < 77,1 and j < 112; 

(ii) all Qifl and Qo,j, where i < ni and j < n 2 \ and 

(iii) the complete energy spectrum of the [n\, 712) system. 
The Q ni .o and Qo.n 2 characterize the single-component 
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Fermi gas and are calculated in Appendix [S] In calculat- 
ing the Q nit n 2 , it is convenient to express the total energy 
fi{ni,n 2 ) ag a sum f ^ ne cen f er f mass energy Ecm and 

the relative energy E^' n2 \ E^ n ^ = E CM + E^ 1 '™ 2 } 
In Eqs. (fT0]) - (fT2)) . the center of mass contribution cancels 
the Qi,q term in the denominator [39|. 

We now evaluate the 6 nii „ 2 at unitarity, i.e., for l/a s — 
0. The calculation of 61,1 is straightforward [39l |. The 
relative energies of the (ni, n 2 ) = (1, 1) system in a har- 
monic trap can be solved exactly using zero-range s-wave 
interactions (48[. The relative energies of the (1,1) sys- 
tem at unitarity are E~^^ = (2q + l/2)hio for I — and 

E^ = {2q + l + 3/2)Huj for I = 1,2,..., where q is a 
non-negative integer. As the relative energy spectrum for 
I > is identical to the single particle energy spectrum, 
we have that Qi.i/Qi.o is equal to Qo,i for states with 
I > 0. The remaining sum, 



bi. 



E 

9=0 



-(2g+l/2)£> 



-(2q+3/2)tD 



(13) 



includes all s-wave states and can be calculated analyti- 
cally [H, 



61,1 = 



D/2 



1 



(14) 



where uj denotes a dimensionless inverse temperature, 
Cj = hui I '(ksT). Equation (fT4|) shows that bi^i is inde- 
pendent of the mass ratio. This is a direct consequence 
of the fact that the relative two-body energy spectrum 
is independent of the mass ratio. Furthermore, like the 
single-component virial coefficients b n (see Appendix [A]) . 
&i 1 is an even function in Hi. 

To calculate b\$ and 62,1, we need the complete relative 
energy spectrum of the (1,2) system (1 light atom and 2 
heavy atoms) and the (2, 1) system (2 light atoms and 1 
heavy atom), respectively. The three-body energies are 
conveniently expressed in terms of the quantity R, 



1/ k for a majority of light species 
k for a majority of heavy species. 



(15) 



At unitarity, the relative three-body energies for inter- 
species s-wave zero-range interaction potentials can be 

written as £ , r ( " 1 1, " 2) = (2q + s hv + l)tw [H], where 
q = 0, 1, . . ., and where the si >v denote non-integer val- 
ues. The si lV depend on the particle symmetry and R, 
and are obtained by solving the five-dimensional hyper- 
angular Schrodinger equation [49|-|5l| . The hyperangular 
eigenvalues were first obtained for three identical bosons 
by Efimov 33], but have since been extended to any par- 
ticle symm etry and mass ratio by both Efimov (34 . |35| 
and others p9M51j . Application of Ref . [51( to the two- 
component system with R, relative angular momentum I, 
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FIG. 1: (Color online) Hyperangular eigenvalues 
function of R for (a) I — 0, (b) I = 1, and (c) I = 2 at 
unitarity. Note that R is shown on a log scale. Solid lines 
show the Si,„ obtained by solving Eq. (|16[) . dashed lines show 
the limiting expressions, Eqs. (|17[) - (|20| ). for k < 1, and thin 
dotted lines show the non- interacting sf*. 



and parity (—1)' at unitarity yields 
V^T(l + 3/2) 



2 ^(l + ^-^,l 



I + 81 v 



(16) 



Z+3 1 

2 '(« + l) 2 



where T is the gamma function and 2-F1 is the hyper- 
geometric function. Equation (|16j) has a spurious root 
at So,;/ = 2 for all re; this root needs to be removed "by 
hand". 

Solid lines in Fig. [1] show the s; )W as a function of R 
for the three lowest relative angular momenta, i.e., for 
I = 0, 1, and 2. The hyperangular quantum number v, 
v = 0, 1, . . ., counts the number of times that the slope 
of si M changes sign. The thin dotted lines indicate the 
non-interacting limits, Sq v = 2v + 4 for I = and sfl = 
2v + I + 2 for I > 0. The dashed lines show an expansion 
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of Eq. (fig for R < 1, 
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3^ 



S ,2j « 4j + 2 + V7IJTT) K, 

V07T 



SQ,2j-l « 4j + 2 - + R, 

V07T 



where j = 1, 2, 



and 



Sl„ « 2;/ + / + 1 + 



(-l)' + T(/ + i/ + l) 

Ar(z + 3/2) r(i/ + i) 



(17) 
(18) 
(19) 

(20) 



for I > 0. For / = 0, the lowest eigenvalue varies quadrat- 
ically with R in the small R regime. The higher lying si >v 
values for I — appear in pairs as R — > 0, with a splitting 
that is linear in R. For I > 0, in contrast, all eigen- 
values are non-degenerate as R — > 0. The lowest I = 1 
eigenvalue equals 2 for R = 0, crosses 1 for k « 8.62, 
and becomes purely imaginary for R ss 13.61. The lat- 
ter point indicates the onset of Efimov physics [28| . For 
8.62 < R < 13.61, non-universal three-body states can 
exist ^28l - l32l | . Similar behavior is found for the lowest 
si >v values for higher odd angular momenta. The si tt , 
values for even I are greater than 1 for all R, indicating 
the absence of non-universal and Efimov physics in the 
even I angular momentum channels. For I = 1,3,5,... 
[see Fig.QTb) for 1 = 1], the si iV with v > begin one unit 
below the non-interacting value in the R — limit and 
decrease to the next lower non-interacting value in the 
limit R — y oo . For I = 2, 4, 6, . . . [see Fig. QTc) for 1 = 2], 
the si yV begin one unit below the non-interacting val- 
ues in the R — limit and approach the non-interacting 
values in the limit R — > oo. In the limit R — > oo, all 
states save those that describe Efimov physics behave 
like non-interacting states. The two masses are so heavy 
compared to the single light particle that even infinitely 
strong interactions cannot mediate a lowering of the en- 
ergy. 

Now that we have the si yV , we can calculate Qi,2 
and Q2,i, and thus b\ 2 and 62,1- We restrict ourselves 
to the regime where Efimov physics is absent (i.e., we 
limit ourselves to R < 13.61) and we assume that the 
three-body system behaves fully universal (i.e., we as- 
sume the absence of three-body resonances in the regime 
8.62 < R < 13.61). Under these assumptions, the virial 
coefficient 62.1, Eq. ([12")) , can be written as 



6 2> l =£>,!(/), 



(21) 



1=0 



where 



mo = £E( 2I+1 ) 

i/=0 q=0 



-(2q+s, iV + l)C 



. e -(2g+s™i+l)« _ e -(2g+l/2+2i/+i+3/2)<D 



+e 



-(2g+3/2+2i/+i+3/2)i> 



(22) 



The first and second terms in the square bracket on the 
right hand side of Eq. (f2"2"j) arise from the Q2,i/<9i,o and 
Q2,o terms, respectively, while the third and fourth terms 
in the square bracket on the right hand side of Eq. (|2"2")l 
arise from the — &1.1Q1.0 term. For / > 0, the second and 
fourth terms cancel, while the third term can be rewritten 
in terms of s^*. Performing the sum over q analytically, 
we find 



M' > °) 



(23) 



(2/ + l)f](e- ( ^ +1 > 



v=0 



For I = 0, the second and fourth terms do not completely 
cancel and contribute, after the sum over q is done an- 
alytically, a term proportional to exp(— 3w). Combining 
the first and third terms in the square bracket on the 
right hand side of Eq. (|2"2"|) , we find 



b 2 ,i(l - 0) 



e 2w - 1 



e -3u, _ e -2u 



(24) 



OO 

u 

u=0 



Since the mass ratio dependence only enters through the 
si fU , Eqs. ([2"T1) - (|2"4"1) remain valid for bi t 2- In practice, we 
calculate the first 10,000 s;.^ values for each I, I — — 50, 
and determine the 61,2 (0 and 62.1 (0 using Eqs. (|2U)) and 
(|24l) . We find that we obtain better convergence if we 
employ the same cutoff on the terms that involve Si v 
and sfl than if we perform the sum over v that involves 
the non-interacting sfl values analytically. 

Figure [5] shows the &2.i(0 for the first five angular mo- 
menta for (a) K = 1 and (b) R = 6.67. The 62,1 (0 are 
negative for even I and positive for odd I for all Co. The 
leading contribution for any mass ratio comes from the 
solution with v = 0. For odd I, we have s;,o + l < s fo [ scc 
Fig. IHb) for I = 1], implying that the 62,1(0 are positive 
for all u and R. For even I, we have s;.o + 1 > sfo I see 
Fig.pTc) for I — 2], implying that the &2,i(0 are negative. 
For I — 0, the — exp(— 2Q) term gives the leading contri- 
bution, implying that &2,i(0 is negative for all Co and R. 
A similar analysis can be conducted for the bi^(l)- The 
alternating sign of the 61,2(0 and the 62,1 (I) implies that 
the full virial coefficients 61,2 and 62,1 converge like an 
alternating series as more I terms are included. We find 
that the convergence rate decreases as R increases. 

Figure [3] shows the virial coefficients 62,1 for R = 2, 3, 4 
and 6.67 (from bottom to top) as a function of Co. The 
virial coefficient 62,1 for R = 2 and 3 is negative for all 
Co. The virial coefficient for R — 4 is negative for large 
Co, changes sign at Co sa 1.67, and is positive for small Co. 
The virial coefficient for R — 6.67 changes sign at a much 
lower temperature, i.e., at Co ss 5.59. The sign change of 
62.1 can be attributed to the fact that the negative I — 
contribution increases slower in magnitude than the pos- 
itive I = 1 contribution with increasing R (see Fig. [5]) , so 
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FIG. 2: (Color online) Angular momentum contributions 
f>2,i(i) for Z = — 4 for (a) R = 1 and (b) R = 6.67 as a 
function of u> at unitarity. Solid, dotted, dashed, dash-dotted, 
and dash-dot-dotted lines show 62,1 (0 for 1—0, 1, 2, 3, and 4, 
respectively. 
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FIG. 4: (Color online) Circles show the universal parts b^{\ 

(R < 1) and &2°i (S > 1) of the virial coefficients as a function 
of k. The inset (the axes use the same labels as the main 
panel) shows the first non-universal corrections and o£ \ 
(squares) and the second non-universal corrections b{ 2 an d 

"21 (triangles) as a function of R. The lines are guides to the 
eye. 




CO 



FIG. 3: (Color online) Solid, dashed, dash-dotted, and dash- 
dot-dotted lines show the virial coefficient 62,1 as a function 
of Co for R = 2, 3, 4 and 6.67. The thin dotted line at is a 
guide to emphasize the sign change of 62,1 for R = 4 around 
u ~ 1.67. 



that the positive / = 1 contribution dominates for suffi- 
ciently large R. In the high temperature limit, b 2 .i first 
becomes positive for It w 3.11. In the low temperature 
regime, 62,1 remains negative for R < 8.62 and changes 
sign when si t o — 1, i.e., the R value beyond which three- 
body resonances can appear [28l432| . We find that the 
virial coefficient 61,2 is negative for all R and G). We note 
that a general framework for calculating 61.2 and 62,1 for 
unequal-mass systems is given in Ref. [52j; however, this 
reference did not provide numerical values for and 
62.1. The fact that the third order virial coefficient 62,1 



changes sign as R increases suggests that the higher or- 
der virial coefficients might also change sign. Through- 
out this paper, we restrict our analysis to systems with 
m + n 2 < 3. 

In the high temperature limit, we Taylor expand the 
b ni ,n 2 (n-i + n 2 = 3) [Hj], 

bn u n 2 = ^ + + bi nL* 4 + (25) 

To our calculated accuracy, we find that the odd powers 
vanish, i.e., the virial coefficients 61,2 and 62,1 appear to 
be even functions in w, similar to their single-component 
counterparts. For the equal-mass system, the expres- 
sion (|25|) was considered in Ref. [39|. The 6n?, na term 
is independent of uj and thus universal. The higher or- 
der corrections 6^™ 2 > °'' an d ^sTi^ are non-universal, since 
they depend on the harmonic trapping potential through 
uj. The supplementary material [53| tabulates the 6^ 
and &| 1, 71 — 0, 2,4 and 6, as a function of R. The b^l 

and &2 1 with n = 0, 2, and 4 are shown in Fig. 2] as a 
function of R, where R is shown on a log scale. The univer- 
sal virial coefficient b 2 °\ has an infinite slope at R = 13.61 
where Efimov physics sets in. For R > 13.61 (not shown), 

6 2 °i is expected to diverge due to the contributions of the 
infinite sequence of Efimov trimers. In the limit R — ¥ 0, 
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the 61,2(0 with I > vanish and we have 

_ p 2ui 

6i, a (K = 0) 



(e Q + l) 2 (e 2Q + 1) 
-ui 



32 



64 



(26) 



III. THERMODYNAMICS OF THE TRAPPED 
MASS-IMBALANCED TWO-COMPONENT 
FERMI GAS 

This section considers the thermodynamics of the har- 
monically trapped two-component Fermi gas at unitarity 
as a function of n. For the equal-mass case, we refer the 
reader to Refs. [39l - l43j . We measure our temperature T 
in units of the semi-classical Fermi temperature Tp of a 
single-component Fermi gas with N/2 particles of mass 
mi, T F = [6(N/2)] 1 / 3 hw/k B Correspondingly, our 

energy is measured in units of Ep, where Ep — kp,Tp. 
Our starting point is the thermodynamic potential fi( 2 \ 
Eq. ([S]), with ni + ri2 < 3 (i.e., including bi 1, b\ 2, and 

&2,l). 

The two coupled number equations, derived from Ni — 
-dCl^/dfii (i=l and 2), are 

= fi(zi) + h :1 z 1 z 2 + 2b 2 ^z 2 z 2 + b ia z\z\ (27) 



Qi,o 



and 



N 2 

Qi,c 



i(z 2 ) + 6i,iziz 2 + 62,1^122 + 26i )2 2iz|. (28) 



Throughout, we use the full temperature dependent ex- 
pressions for the virial coefficients, and 



/i( z i) = X! n ^ifi z i 



m=o 



and 



/a (^2) = X! n2 &0,n 2 Z : 



"2 

2 • 



(29) 



(30) 



n 2 =0 



Equations ([29]) and ([30]) can be employed if z, ; < e 3i / 2 , 
i = 1,2. For larger z i7 useful parameterizations of 
Eqs. O and §U§ are 



fi{Zi) ~ -Li 3 (-Zi) 

+ \u 2 [-Ln(l + Zi) - U 3 (- Zi )] 



1 ,-.4 
1920^ 



17 z t 



(31) 

30Ln(l + z i )-13Li 3 (-z i 



where Ln is the natural logarithm function and Li is the 
polylogarithm function. For the temperatures considered 




FIG. 5: (Color online) Trapped system at unitarity. Solid and 
dashed lines (dash-dotted and dotted lines) show the fugaci- 
ties z\ and Z2, respectively, for « = 1 [k = 13) as a function 
of the polarization P for TV = 5 x 10 7 and (a) T/Tf = 1 and 
(b) T/T F = 1/2. 



in this paper, terms of order uj e and higher are negligi- 
ble. Equations ([2"9"]) and (13TJ1) treat the non-interacting 
"reference pieces" fi(zi) and f 2 (z 2 ) as an infinite series 
in the fugacities and the interacting pieces at third or- 
der in the fugacities. Alternatively, one may consider 
truncating the non-interacting and interacting pieces at 
the same order in the fugacities. We have opted for the 
former approach for two reasons: (i) In the limit that 
P = 1 or —1, we recover the exact results of the single- 
component Fermi gas. (ii) In the regime where the z% 
are small, i.e., in the regime where the virial equation 
of state up to third order is expected to provide reliable 
results, we find fairly small differences [e.g., < 0.06% in 
Fig. [5{a)] between the approaches that treat the non- 
interacting piece at infinite order and third order in the 
fugacities. For a given temperature T, polarization P, 
and total number of particles N, we solve Eqs. ([2"T]) and 
([2"8"]l self-consistently for z\ and z 2 . 

Figure [5] shows the fugacities z\ and z 2 at unitarity 
as a function of the polarization for (a) T /Tp = 1 and 
(b) T/Tp = 1/2 for k — 1 (solid and dashed lines) and 
K = 13 (dash-dotted and dotted lines). The fugacities for 
T/Tp = 1 are smaller than 0.35 for both mass ratios con- 
sidered, suggesting that the virial equation of state can 
likely be trusted for T/T F > 1. For T/T F = 1/2, in con- 
trast, the fugacities are as large as 3.5, suggesting that 
the description is at best qualitatively correct. On the 
scale shown, the fugacities for T/Tp = 1 and k — 1 are 
nearly indistinguishable from those at the same tempera- 
ture but higher mass ratio. This makes sense intuitively, 
since the fugacities are small in the high temperature 
regime, thereby suppressing higher order terms in the 
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FIG. 6: (Color online) Free energy per particle F/N (N = 
5 x 10 7 ) as a function of the polarization P for (a) T/Tf — 1 
and (b) T/Tf = 1/2 for the trapped system at unitarity. 
Alternating solid and dotted curves, from top to bottom, show 
the free energy for n — 1 — 13 in units of 1. The dashed lines 
connect the (F eq ,P eq ) points for different re. 



virial expansion. At lower temperatures [see Fig. EJb)], 
the higher order virial coefficients become more impor- 
tant, as evidenced by the fact that the fugacities show a 
notable dependence on the mass ratio. 

In the following, we use the virial equation of state to 
determine the free energy F, F = flW + fiiNi + ij, 2 N 2 , 
the entropy S, S = — dfl&'/dT (calculated for fixed w, 
Hi and and the energy U, U — F + TS, as a function 
of the mass ratio re, the polarization P, and the tempera- 
ture T at unitarity. Solid and dotted lines in Fig. [5] show 
the free energy per particle F/N as a function of the po- 
larization for (a) T/Tp = 1 and (b) T/T F = 1/2. In both 
panels, the alternating solid and dotted lines, from top to 
bottom, correspond to mass ratios re = 1, 2, . . . , 13. For 
the equal mass system, the minimum of the free energy 
occurs for P = 0, i.e., for equal numbers of spin- up and 
spin-down fcrmions. As re increases, the minimum of the 
free energy moves to positive polarizations, i.e., to a ma- 
jority of heavy particles. We refer to the (F, P) values at 
which the free energy is minimized as (F eq ,P eq ). Dashed 
lines in Figs. [B^a) and [B^b) connect the (F eq ,P eq ) val- 
ues for different re. For both temperatures, P eq increases 
with increasing re, which can be understood by realiz- 
ing that the three-body system consisting of two heavy 
atoms and one light atom has a lower energy than the 
three-body system consisting of one heavy atom and two 
light atoms. Moreover, for a given re, P eq increases with 
decreasing temperature. This trend makes sense, since 



we expect interaction effects to become more important 
as the temperature decreases. 

To obtain an analytical expression for P eq , we take 
the derivative of F with respect to P at fixed temper- 
ature and fixed number of particles. We find that the 
extremum of P occurs when the two fugacities are equal 



(i.e., zi 



Z-2 



j). This restriction allows for a 



straightfoward calculation of z eq , namely 



N 



Q 



1.0 



fl(Zeq) + fiiZeq) 



3(62,1 + bl :2 )ze q 

(32) 



can be solved self-consistently for fixed N and T. Up to 
third order in z eq , P eq can be written as 



Pr, 



N 



/i,o 



(6 2 ,i - b la )z. 



(33) 



P eq is equal to zero for all re if the virial equation of state 
terminates at the second order in the fugacity. At the 
third order, we find P eq = for re = 1, since &i ; 2 equals 
62,1 in this case. For re > 1, however, P eq is positive 
since 62,1 > 61,2 for all re > 1. Interestingly, the sign 
change of 62,1 in the high temperature regime for re w 
3.11 does not lead to a discontinuity or a sign change of 
P eq since is finite and negative for all re. For small 
deviations from re = 1, we find that P eq and F eq change 
linearly and quadratically, respectively, as a function of 
re for fixed temperature and fixed number of particles. 
We find that d 2 F/dP 2 , calculated for fixed T and N, 
is greater than zero for all polarizations and mass ratios 
considered, which implies that the system is stable at this 
level of approximation. 

To discuss the dependence of the entropy S and the 
energy U on the temperature and the polarization, we 
focus on the mass ratio re = 6.67, i.e., on the experimen- 
tally realizable K-Li system at unitarity |2fj| - l2a |. Thin 
dotted and solid lines in Fig. Eta) show the free energy 
per particle as a function of the polarization for T/Tp = 
1, 1.2, ...,1.8 (from top to bottom). The dashed line con- 
nects the (F eq ,P eq ) values for different T (but fixed re 
and N). The vertical solid line marks the value P eq for 
T/T F = 1, where F/N = -2.82E F , S/N = 5.73k B , and 
U/N = 2.QIE F . Thin solid and dotted lines in Fig. [TJb) 
show the entropy per particle S/N as a function of the po- 
larization for T/T F = 0.95, 0.96, . . . , 1.04 (from bottom 
at P to top at P « 0.5). Similarly, thin solid and dot- 
ted lines in Fig.JTJc) show the energy per particle U/N as 
a function of the polarization for T/Tp — 0.99, 1, . . . , 1.06 
(from bottom at P ~ 0.5 to top at P w 0.1). The dashed 
lines in Figs. Elb) and [7jc) connect points of constant 
U/N = 2ME F and S/N = 5.73k B , respectively (i.e., 
the values of U/N and S/N at P eq for T/Tp = 1). As 
expected, the entropy for a fixed energy and the energy 
for a fixed entropy are concave and convex functions of 
the polarization, respectively. It can be seen that the 
entropy per particle S/N is maximized at P eq while the 
energy per particle U/N is minimized at P eq . The fact 
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FIG. 7: (Color online) Thermodynamic observables F/N, 
S/N, and U /N as a function of the polarization P for k = 6.67 
and N = 5 x 10 7 for the trapped system at unitarity. Alter- 
nating dotted and solid lines show (a) the free energy per par- 
ticle F/N for T/T F = 1, 1.2, . . . , 1.8 (from top to bottom), (b) 
the entropy per particle S/N for T/T F = 0.95, 0.96, . . . , 1.04 
(from bottom at P ~ to top at P w 0.5), and (c) the energy 
per particle U/N for T/T F = 0.99, 1, . . . , 1.06 (from bottom 
at P ~ 0.5 to top at P ~ 0.1). In panel (a), the dashed 
line connects (F eq ,P eq ) values for different T. In panels (b) 
and (c), the dashed lines show S/N for U/N = 2.91.Ef and 
U/N for S/N = 5.73fcs, respectively. In all three panels, the 
vertical lines indicate the value of P eq for T/T F = 1. 



that —d 2 S/dP 2 (d 2 U/dP 2 ) is greater than zero for fixed 
U and N (for fixed 5 and N) is consistent with the fact 
that the system is stable as concluded earlier from the 
fact that d 2 F/dP 2 is greater than zero for fixed T and 
N . Lastly, we note that the maximum of the isotherm 
S/N moves to smaller P with decreasing T and that the 
minimum of the isotherm U/N moves to larger P with 
decreasing T. 

The thermodynamic quantities shown in Figs. [S][7J have 
been obtained for TV = 5 x 10 7 particles. For fixed T/Tp, 
the actual temperature T decreases with decreasing N, 




FIG. 8: (Color online) Percent difference between P eq for iV = 
5 x 10 7 and P eq for iV = 500 as a function of the mass ratio k 
for the trapped system at unitarity. Solid, dashed and dotted 
lines show the percent difference for T/T F — 1/2, 1 and 2, 
respectively. 



since Tp decreases with decreasing N. Correspondingly, 
u! increases with decreasing N (again, assuming fixed 
T /Tp). Since the leading order non-universal corrections 
scale with Co 2 [see Eqs. (f^5|) and (jBTjl ] . the non-universal 
corrections are expected to be negligible for large N, but 
not for very small N. Figure [5] shows the percentage dif- 
ference between P eq for N — 5 x 10 7 and P eq for N = 500 
as a function of k for T/Tp = 2, 1 and 1/2. The non- 
universal corrections are largest, i.e., on the order of 2%, 
for k = 1 and T/T F = 1/2 [solid line in Fig. 0. For 
larger N, the non-universal corrections are even smaller, 
implying that Figs. [5J7J show essentially universal, TV- 
independent thermodynamic quantities of the trapped 
two-component Fermi system with equal trapping fre- 
quencies. 

We conclude this section by noting that P eq is quite 
small for the K-Li system down to T/Tp — 1, but then 
changes more rapidly as the temperature is lowered fur- 
ther. The dependence of P eq on the mass ratio k may be 
enhanced by switching to a trap geometry for which the 
single-particle energies are dependent on the mass mi or 
TO2- One such system is, as detailed in the next section, 
the homogeneous system, which favors much larger polar- 
izations than the trapped system for comparable T/Tp. 



IV. THERMODYNAMICS OF THE UNIVERSAL 
HOMOGENEOUS MASS-IMBALANCED 
TWO-COMPONENT FERMI GAS 

A. Local Density Approximation 

This section shows how the homogeneous system can 
be related to the harmonically trapped system via the 
local density approximation [33 . |55| . In the high tern- 
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(2) 

perature limit, the thermodynamic potential of the 

homogeneous system can be written as 



rt 2) = nft 

horn l.hom 



^2, horn + ^^hom, 



(34) 



where 



n 



(i) 

l.hom 



(i) 

2, horn 



v 

A «ii m=i 
V °° 



z?\ (35) 



om "\ i 



(36) 



om ^2 ' 



and 



Afihom — kpjT—^ S ' ^ ' 6 r ii ,n 2 ,hom z i 1 ^2 2 ■ (^7) 
V m=liij=l 

Here, the & nil n 2 ,hom denote the virial coefhcients of the 
homogeneous system, and 21 and zi are the fugaci- 
ties of the homogeneous system. The virial coefficients 
frni,n 2 ,hom are defined by Eqs. (fTUj) - (fl2")) with the Q ni ,n 2 
interpreted as the partition functions of the homogeneous 
systems. Equations (f3"4" |) - ([3"T)) are analogous to the high 
temperature limits of Eqs. dU)-©. In Eqs. ([M ll -lpTF ]) . we 
used that the high temperature limit of <3i,o,hom is given 
by V/X 3 ni [IBj], where V is the volume of the system and 
where A mi is the thermal de Broglie wavelength of a par- 
ticle of the heavy species, 



00 00 



Ami 



' 2irh 2 
miksT 



(38) 



In the local density approximation |39|, [55| , the ther- 
modynamic potential fl^ of the trapped system is given 



by a weighted average of fi[ 2 ' m (r) over all space, 



Jd 3 r 



(39) 



i.e., we assume that the thermodynamic potential at 
each location r in the trap can be approximated by 
the corresponding bulk value. This implies that the 
Zi in Eqs. ([3"5 j) -([3"T f are replaced by Zj(r), Zi{r) — 
exp[jj,i t i oc (r) I (feflT)], where the local chemical potentials 
(f) are defined as 



Hi,ioc(r) = ^ 



1 2 -Q, 

^rriiUJ r . 



(40) 



The integrations on the right hand side of Eq. (|3"9")l 
are straightforward. We first perform the integration in- 
volving Af2hom(^*)- The integration in the denominator 
introduces a factor of V, which cancels with the factor 
of V in the numerator. This leaves, after performing the 
angular integrals, 



Afi = -47r^ £ E Ws-wr^x (41) 

™i ni=lna=l 

poo 

/ exp [— i(mini + rri2n2)uJ 2 r 2 / (kBT)]r 2 dr . 



00 00 



Performing the remaining integration and expressing the 
result in terms of the mass ratio k, we find 



-k B T 

OO OO 



k B T 



%i ,n 2 ,hom 



1 1 K + n 2 /nf/ 2 1 

rii- 1 no — 1 



(42) 



The prefactor [k bT / '(hu)} 3 coincides with the high T 
limit of Qifl for the trapped system. Comparison of 
Eqs. (|9|) and (j42|) allows us to relate the universal virial 
coefficients bn},n 2 of the trapped system to the virial co- 
efficients 6 nil n 2 .hom of the homogeneous system, 



&n 1 ,n a ,horo = ("1 + n%/Kf /2 b^ n2 



(43) 



Considering the integrations over ^i\ wu Sf) an d 



^ihomW' one nn ds analogous relationships between 
&rn,o,hom and o^l fi , and between & ,n 2 ,hom and 6^ 2 . 



B. Universal thermodynamics of the homogeneous 
system 

This section considers the universal, iV-independcnt 
thermodynamics of the homogeneous two-component 
Fermi gas at unitarity as a function of n. We 
measure our temperature T in units of the semi- 
classical Fermi temperature Tp of a single-component 
Fermi gas with N/2 particles of mass mi, fcsTp = 
[6Tr 2 (N/2)/V] 2 / 3 h 2 /(2m 1 ) [H. Notice that the Fermi 
temperature scales inversely with the mass mi. This 
implies that the Fermi temperature of N/2 light parti- 
cles with mass TO2 is n times larger than the Fermi tem- 
perature of N/2 heavy particles. The different Fermi 
temperatures of the single-component heavy and light 
particles are, of course, a direct consequence of the mis- 
matching Fermi surfaces of unequal-mass Fermi gases or, 
equivalently, a direct consequence of the fact that the 
single-particle energies of a mass mi particle in a box 
with periodic boundary conditions scale with l/m*. The 
high temperature virial equation of state of the homo- 
geneous system can be applied down to T/Tp about 1 
for k = 1. For unequal masses, the applicability regime, 
using the heavy mass to define Tp (see above) , is limited 
to T /Tp > k. 

The coupled number equations of the homogeneous 
system, up to third order, read 

= - Li 3/2 (-Zi) + b 1A 



and 



V 



3/2-^2 A mi t- ( \ 1 3/2l 

K 1 y = -Ll 3/2 (-Z2) + K ' 0i,i,hom«l«2 

+ K 3/2 &2,l,hom2^2 + 2n 3/2 b 1 ,2,ho™ z l z 2- ( 4 5) 



hom^l^ 



2&2,l,hom2l 22 + bl^horo^l^ ( 44 ) 
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FIG. 9: (Color online) Fugacities z\ and 22 of the homo- 
geneous system as a function of the polarization P for (a) 
k = 1 and T/T F = 6.67, (b) k = 6.67 and T/T F = 6.67, (c) 
k = 1 and T/T F = 2, and (d) ft = 2 and T/T F = 2. Dotted 
and dash-dotted lines show z\ and 22, respectively, for the 
unitary system obtained by solving Eqs. (|44p and (|45p self- 
consistently. For comparison, solid and dashed lines show z\ 
and 22, respectively, for the non-interacting system. 



Dotted and dash-dotted lines in Fig.Oshow z\ and Z2, re- 
spectively, as a function of the polarization P for the uni- 
tary system obtained by solving Eqs. (|4"4")) and (f4"5j) self- 
consistent ly for fixed n and T. The panels show the fu- 
gacities for (a) K = 1 and T/Tp = 6.67, (b) k = 6.67 and 
T/T F = 6.67, (c) K = 1 and T/T F = 2, and (d) « = 2 and 
T/T F = 2. For the equal mass case, Eqs. (gt]) and (gSJ) 
are symmetric and thus z\ and Z2 are symmetric about 
P = for any temperature [see Figs.^a) and^c)]. For 
unequal masses, in contrast, Eqs. (|44l) and (|45|) are not 
symmetric and the crossing of the fugacities Z\ and z^ is 
shifted to positive polarizations [see Figs.|Hfb) and^d)]. 
For comparison, solid and dashed lines in Fig. |H] show the 
fugacities Z\ and z%, respectively, for the corresponding 
non-interacting systems. Even for the non-interacting 
system, the fugacities z\ and z-i are not symmetric with 
respect to P = 0; this is a direct consequence of the k 3 / 2 
factor in Eq. (|45|) . For — 1 < P < 1, the fugacities for 
the non-interacting systems lie above the respective fu- 
gacities for the unitary systems. For k > 1, this implies 
that the unitary system has a smaller equilibrium polar- 
ization P eq than the non-interacting system for the same 
T and k. Intuitively, this makes sense since we expect 
that the attractive interactions, at least for mass ratios 
not too much larger than 1, push the system toward an 
equal mixture of heavy and light particles. 

Figure [TO] shows the free energy per particle F/N as 
a function of the polarization P for TjTp = 10 and 
three different mass ratios k. The solid lines are cal- 
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FIG. 10: (Color online) Free energy per particle F/N of the 
homogeneous system as a function of the polarization P for 
T/Tf = 10 and mass ratios k = 1, 2, and 6.67. The curves are 
grouped in sets of two, each labeled by their corresponding 
mass ratio. The solid lines are calculated using the virial 
equation of state at unitarity. For comparison, the dotted 
lines show F/N for the non-interacting system. 



culated using the virial equation of state of the unitary 
system. For comparison, dotted lines show F/N for the 
non-interacting system. Since the free energy is a con- 
cave function with respect to P for fixed N and T, the 
equilibrium polarizations are stable minima. 

To better understand the interplay of the second and 
third order virial coefficients for the homogeneous sys- 
tem, we calculate z eq analogously to Eq. (|32[) and find 



P eq for the homogeneous system to be 



P NX ^ 

r eq — 



V 



- Ll3/2(-Zeg) 



K 3/2 



[b 



2,1. horn 



2, ho 



(46) 



Eliminating the polylogarithm function, P eq can be 
rewritten as 



P,, 



t 3 / 2 - 1 
.3/2 + 1 



l.hom 



k 3 / 2 + 1 N\z mi Zeq ~ 



1 - 2k 3 / 2 



(n) 



I ' 2,1 



,3/2 



horn 
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(47) 



2V 3 
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As expected, the z eq term depends only on the mass ra- 
tio k and is greater than zero for k > 1; in the limit 
k — y 00, the non- interacting system prefers to have only 
heavy particles, i.e., P eq — > 1. The z 2 term of the uni- 



eq 

tary system, in contrast, is negative for all k > 1. This 
implies that the second order virial coefficient of the uni- 
tary system acts to decrease P eq compared to the P eq of 
the non-interacting system. Interestingly, the z 3 q term of 
the unitary system is positive for k < 5.40 and negative 
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FIG. 11: (Color online) Equilibrium polarization P eq of the 
homogeneous system at unitarity as a function of the mass 
ratio k for T/Tf = 10. The dashed and solid lines show 
P eq calculated using the virial equation of state up to second 
and third order, respectively. The dashed and solid lines are 
nearly indistinguishable on the scale shown. For comparison, 
the dotted line shows P eq for the non-interacting system. The 
solid line in the inset shows the difference between the P eq s 
calculated using the virial equation of state up to the third 
and second orders (the dashed line shows zero as a reference). 

for k > 5.40, independent of temperature. This implies 
that the third order virial coefficients act to increase (de- 
crease) P eq for 1 < k < 5.40 (k > 5.40) compared to 
the P eq calculated using the virial equation of state up 
to second order. 

Figure QT] shows the equilibrium polarization P eq as a 
function of n for T/Tp = 10. The dashed and solid lines 
show P eq calculated using the virial equation of state of 
the unitary system up to second and third order, respec- 
tively. The second and third orders are nearly indistin- 
guishable on the scale shown. For comparison, the dotted 
line shows P eq for the non-interacting system. The inset 
of Fig. [TT] shows the difference of P eq calculated up to 
the second order subtracted from P eq calculated up to 
the third order. Note that the zero crossing occurs at 
k = 5.56, and not at k — 5.40, since the equilibrium fu- 
gacity determined from the virial equation of state up to 
third order is slightly smaller than that determined from 
the virial equation of state up to second order. 

Figure [12] shows the equilibrium polarization P eq as a 
function of T for k = 6.67. The dashed and solid lines 
show P eq at unitarity calculated using the virial equation 
of state up to second and third order, respectively. The 
second and third orders are nearly indistinguishable on 
the scale shown. For comparison, the dotted line shows 
P eq for the non-interacting system; it is independent of 
T and equals 0.89. The inset of Fig. shows the dif- 
ference of P eq at unitarity calculated up to the second 
order subtracted from P eq at unitarity calculated up to 
the third order. Since the mass ratio used is larger than 
5.40, the third order correction acts to further decrease 



1 1 — i — i — i — i — | — i — i — i — i — | — i — i — i — i — | — i — i — i — r 




FIG. 12: (Color online) Equilibrium polarization P eq of the 
homogeneous system at unitarity as a function of the tem- 
perature T for ft = 6.67. The dashed and solid lines show 
P eq calculated using the virial equation of state up to second 
and third order, respectively. The dashed and solid lines are 
nearly indistinguishable on the scale shown. For comparison, 
the dotted line shows P eq for the non-interacting system. The 
solid line in the inset shows the difference between the P eq 's 
calculated using the virial equation of state up to the third 
and second orders. 

P eq for T/Tp > 5. Note that the virial equation of state 
is likely unreliable for T/Tp < 6.67. As the tempera- 
ture decreases, interactions become more important and 
the deviations between P eq for the unitary and the non- 
interacting system increase. 

V. CONCLUSIONS 

This paper considers the high temperature virial equa- 
tion of state of unequal-mass two-component Fermi gases 
at unitarity. For the trapped system with equal frequen- 
cies and K > 1, we found that the equilibrium polariza- 
tion P eq at unitarity is relatively small for the temper- 
atures considered. The small deviations from P eq = 
are introduced by the third order virial coefficients. For 
the homogeneous system, in contrast, we found compar- 
itively large P eq for k > 1. The value of P eq is de- 
termined by an intricate interplay between the second 
and third order virial coefficients. For the cases inves- 
tigated, P eq increases with decreasing temperature for 
the trapped system while P eq decreases with tempera- 
ture for the homogeneous system. The fourth-order virial 
coefficients 63^1 and 62,2 for unequal-mass systems are 
presently unknown, preventing us to estimate the impor- 
tance of higher-order terms in the virial equation of state. 

At the temperatures considered, the curvature of the 
free energy is positive for all polarizations, thus implying 
the absence of any first-order phase transition. If the 
unequal-mass system is prepared with a polarization P 
different from P eq , then the free energy lies above the 
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value it would have for P eq . This suggests that the system 
prepared with P ^ P eq may phase separate when the 
temperature falls below a certain critical value. 

In future work, it will be interesting to investigate the 
high temperature thermodynamics of unequal mass sys- 
tems when the light and heavy species are confined by 
harmonic potentials with different frequencies via the lo- 
cal density approximation. Moreover, the treatment of 
cylindrically symmetric traps will be relevant for ongo- 
ing experiments. Other future extensions of our work in- 
clude the treatment of unequal mass systems away from 
unitarity and at lower temperatures. 
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FIG. 13: (Color online) Virial coefficients 6„ of the single- 
component harmonically trapped Fermi gas as a function of 
uj. Solid, dotted, dashed, dash-dotted and dash-dot-dotted 
lines show 62,63,64,65 and be, respectively. 



Appendix A: Thermodynamics of the trapped 
single-component fermi gas 

This appendix studies the virial expansion of the 
single-component Fermi gas in a harmonic trap. The 
results for the single-component Fermi gas serve as a ref- 
erence for the trapped s-wave interacting two-component 
Fermi gas. We first derive compact expressions for the 
virial coefficients, and then discuss convergence proper- 
ties of the virial expansion. 

Equation (J7J of the main text defines the expansion 
coefficients b n \ throughout this appendix, we replace n\ 
by n and use a single subscript n to denote the virial 
coefficients b n and the canonical partition functions Q n . 
Partial derivatives of the grand canonical potential fiW 
determine the various thermodynamic quantities. The 
number of particles N of the single-component Fermi gas, 
e.g., can be expressed as a derivative with respect to the 
chemical potential fi, 



N 



= Qi nb ri 



(Al) 



To determine the coefficients b n , we use that the popula- 
tion N(Ej) of a state with a given single-particle energy 
Ej is given by the Fermi-Dirac distribution, 



N(E 3 ) 



z {E 3 -^)/(k B T) + ]_> 



(A2) 



with the constraint that the total number of particles N 
is given by 



n = J2n(Ej). 



(A3) 



The single-particle energies of the three-dimensional har- 
monic oscillator are given by Ej = (j + l/2)hui, j = 
1,2,3..., and have a degeneracy of j(j + l)/2. Us- 
ing the single-particle energies and their degeneracies in 
Eqs. (|A2j) and (|A3j) . we find a relationship between N, 
jj,, and T, 



00 



i(i + i) 



(A4) 



To compare Eqs. (|A4|) and (jAlj) . we Taylor expand Eq. 
4| for small z, z = exp[/i/(fcsT)], 



3=1 



j(j + 1) 



n=l 



(-1) 



n+1. 



(A5) 



The sum over j can be done analytically. Pulling out a 
factor of Qi, we obtain 



N 



e'2 



(n-l)C 



1 



1 



(A6) 



Comparing Eqs. (|A6I) and (|A1I) . the expansion coeffi- 
cients b n can be read off, 



(-1) 



n+1 



e2 



(n-l)u 



1 



- 1 



(A7) 



Figure fTS] shows the first five expansion coefficients as 
a function of uj. In the low temperature limit, the b n 
decay exponentially as exp[— 3(n — l)cD/2]. The sign of 
b n is given by (— l) n+1 for all to and the b n monotonically 
decrease in magnitude as Q increases. The coefficients are 
even functions in cD, i.e., b n remains unchanged as Co is 
replaced by — ui. Expanding Eq. (| A7|) for large T (small 
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FIG. 14: (Color online) Contour plot of the percent difference 
100(z+2&22 2 +3632 3 -A7Qi)/(A7Qi). The dashed line shows 
the percent difference of 1. The solid lines below the dashed 
line show the percent difference in steps of powers of 10 down 
to a value of 0.001, while solid lines above the dashed line 
show the percent difference in steps of 1 up to a value of 7. 
We consider the region where the percent difference is less 
than 1 to be well converged. 



ui), we find [56| 
1 



-1 



+ y^C 17 " 4 _ 30n2 + 13 ) lIj4 + °( q6 ) 



(A8) 



In the limit u> —¥ we have "universal" physics in the 
sense that the virial coefficients b n approach constants, 



b n -> b [ „ [>} as u -> 0, where b ( n> = (-l)" +1 /n 4 [3|. 

Figure PHI shows a contour plot of the percent difference 
between z + 2b 2 z 2 + 3b 3 z 3 and N/Qi [see Eq. ([XI])]. The 
deviation between the exact expressions and the virial 
expansion is smaller than about 1% for z < 1 and is 



below 5% for z w 1.5 for the range of temperatures 
shown. Thus, the virial expansion up to the third order 
provides a qualitatively correct description of the high- 
temperature physics for z < 1.5. Even though we ex- 
panded around small z, Eq. (|A1|) is an analytic function 
of z that can be continued along the positive real axis 
for all z > 1 57]. For the two-component Fermi gas 
with pairwise interactions, in contrast, the convergence 
radius of the virial equation of state, to the best of our 
knowledge, is not known. 

Lastly, one can use Eq. (| A7|) to calculate the Q n for 
n > 1. A compact expression for the Q n , n > 1, can be 
found in the literature [H, H3, HI] , 



Or. 



En 

{m} i=l 



1 



(QihY 



(A9) 



where {m} represents all possible sets of n non-negative 
integers {mi, . . . , m n } that fullfill the constraint 



»=i 



(A10) 



Each set describes a unique way to divide a system of n 
particles into smaller clusters of rrii particles. For n = 3, 
e.g., we have the sets {3, 0, 0}, {1, 1, 0}, and {0, 0, 1}, i.e., 
the three-body system can be thought of as consisting of 
three monomers, of one monomer and one dimer, or of 
one trimer. For completeness, we report the expressions 
for the first few Q n , 



oo oo 



Qi=£5>Z + l)e- 

n=0 1=0 

Qa = hQi + §<9?, 



(2ri+i+3/2)i 



(e" - If 



Q 3 

and 



b 3 Qi+b 2 Ql + ±Ql 



Qa = b 4 Qi + b 3 Qf + \b\Q\ + \b 2 Q\ + ±Qf, 
where we have used that b\ = 1 [see Eq. 



(All) 

(A12) 
(A13) 



(A14) 
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This material tabulates the expansion coefficients 6^ an( i ^2 ™l ( n = 0>2, 4, and 6) of the two-component unequal- 
mass Fermi gas in a harmonic trapping potential at unitarity. The calculation of the virial coefficients 61,2 and 
&2.1 requires the eigenvalues si v [see Eqs. (21)-(24) of the main text]. We calculate the first 10,000 si. u for each I, 
1 = — 50, and each mass ratio considered. To obtain the accuracy reported in Tables U and UH we must maintain 
a higher working precision than machine precision. We use Mathcmatica to determine the si tV with an accuracy of 
about one part in 10 80 (most likely, a notably reduced precision would work just as well). To obtain the b^l and 

&2™i) we Taylor expand Eqs. (23) and (24) of the main text to high order at a finite but small u. The cu is chosen 
such that b± t 2 and 62,1 are converged at that value; note, since we use a cutoff in v and I, we cannot Taylor expand 
around Co = 0. To determine the expansion coefficients for Q — 0, we rewrite the Taylor expanded series in the form 
of Eq. (25) of the main text. For all values listed in Tables HI and ITT1 the error is indicated by the value in parentheses. 
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TABLE I: High temperature expansion coefficients of\, o^\, o^\, 
various mass ratios k. 
67907169/49906069. 



and bj 6 2 for one heavy fermion and two light fermions for 
The odd order coefficients are zero to the calculated accuracy. For the mass ratio of 13.607, we use 
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TABLE II: High temperature expansion coefficients 6 2 i> b% \ , b% \ , and 621 for two heavy fermions and one light fermion for 
various mass ratios k. The odd order coefficients are zero to the calculated accuracy. For the mass ratio of 13.607, we use 
67907169/49906069. 
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